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Abstract 

This paper introduces the fundamental continuum theory governing momentum transport in 
isotropic nanofluidic flows. The theory is an extension to the classical Navier-Stokes equation, 
which includes coupling between translational and rotational degrees of freedom, as well as non¬ 
local response functions that incorporates spatial correlations. The continuum theory is compared 
with molecular dynamics simulation data for both relaxation processes and fluid flows showing 
excellent agreement on the nanometer length scale. We also present practical tools to estimate 
when the extended theory should be used. It is shown that in the wall-fluid region the fluid 
molecules align with the wall and in this region the isotropic model may fail and a full anisotropic 
description is necessary in order to describe this region. 
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I. INTRODUCTION 


Nanoscale devices can now be fabricated with channels where the smallest dimension is 
just a few nanometers [1], and the development of nanofluidic theory [1-3] is more relevant 
than ever. Consider the following example. Perrson et al. [4] used a series of rectangular 
nanochannels with widths ranging from 14 to 300 nm to connect two micro-scale chambers. 
By means of capillary hlling, fluid from one chamber fills up the channels and thus connects 
the two chambers. The hlling rate can be measured for different channel widths and for both 
milli-Q water (hitrated de-ionized water) and an electrolyte solution of sodium chloride. The 
rate did not follow the Washburn equation for channel widths smaller than 100 nm. The 
Washburn equation is based on the classical continuum picture [3] using Poiseuille law of 
huid motion which includes the Newtonian (or macroscopic) shear viscosity. For widths 
larger than 100 nm the Washburn equation correctly predicts the hlling rate. This is in 
accordance with the common understanding that the discrete nature of the huid at small 
scales destroys the continuum picture [2, 6]. In fact, many researchers categorize continuum 
physics as physics on the macroscopic scale, see for example Ref. 7. Several questions 
immediately arise: When exactly does the continuum picture fail? How is this breakdown 
manifested? Does the length scale of the breakdown depend on the specihc problem? Can 
one improve the continuum description such that it applies on small scales? 

Demanding sufficient smoothness of the macroscopic quantities with respect to time and 
position and using a simple statistical argument, Lautrup [7] estimates that the smallest 
volume accessible to the continuum description must contain at least 10^ molecules. This 
corresponds to a length scale of 8-80 nm, depending on the density. For steady hows the 
temporal huctuations can be averaged out and the accessible volume is much smaller. This 
is illustrated in Fig. 1, where we have performed an atomistic simulation (data given by 
blue hlled circles) of a methane huid conhned between two graphene sheets undergoing a 
Poiseuille how. The slit-pore has a width of approximately 3.3 nm and the how is driven 
by an external force held. In this paper the simulations are carried out using the seplib 
library [37]. The classical continuum prediction is plotted as two red lines illustrating the 
maximum and minimum prohles allowed within statistical uncertainty on the Newtonian 
shear viscosity [18]. Only the huid slip velocity at the wall surface is used as a htting 
parameter. For this system the continuum theory gives a satisfactory description of the 
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FIG. 1. (Color online) Comparison between atomistic simulations (blue filled circles) and the 
continuum prediction (red lines) for a methane fluid undergoing a Poiseuille flow. The flow is 
generated by an external force field with magnitude F = 50 TN, pointing along the x-direction. 
The Navier-Stokes equation predicts a velocity profile Ux{z) = pF/{2rjQ){h? — z^) + where 
p = 270 kgm“^ is the mass density, r/o = 9.3 ±0.6 Pa-s the Newtonian shear viscosity, and Uyj = 62 
ms“^ is the fluid slip velocity at the wall surface. The two lines represent the interval associated 
with the standard error in the viscosity [18]. The width of the slit pore is approximately 10 
molecular diameters or 3.3 nm. 


fluid average velocity a length scale of a few nanometers. Apparently, even on these small 
length scales the molecular structure and degrees of freedom can be coarsened into simple 
transport coefficients like the viscosity. For water undergoing a steady flow it has been 
shown by atomistic simulations that the continuum description holds for channel widths of 
just 6-10 nm [1, 9]. These results contrast earlier assumptions about the validity of the 
continuum picture, and the statement that continuum physics is physics on the macroscopic 
scale [2, 7]. Interestingly, it was later argued by Thamdrup et al. [10] that the disagreement 
between the experiment by Persson et al. [4] and the Washburn prediction is due to pinned 
micro-bubles resulting in an increase in hydraulic resistance. 

At some point the classical continuum description will of course break down. To mention 
two examples, Travis et al. [11] showed that for atomic fluidic systems the velocity prohle 
features modulations for conhnements in the order of 5 atomic diameters. Decheverry and 
Bocquet [12] analyzed the effect of thermal fluctuations on mass transport of fluid through 
a nanotube. Interestingly, when the classical continuum theory fails, the dynamics is fre- 
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quently quantified by different transport coefficients compared to those of the bulk system 
and effective transport coefficients are introduced into the continuum constitutive relations 
[13, 14], 

The main point of this paper is that the observation of a breakdown need not be a fail¬ 
ure of the continuum picture itself, but a result of inadequate modeling wherein important 
dynamical processes are not accounted for by classical theories. A very well understood 
example is the effect of the Debye screening layer in electrolyte micro-flows [3]. Two other 
physical mechanisms that become important on the nanoscale are often ignored in the lit¬ 
erature, and this paper will treat these in detail; 

(i) In classical hydrodynamics the fluid’s local rotation is determined uniquely by the 
fluid streaming velocity. One can quantify the rotation from the local angular velocity held 
which is one half the vorticity, that is, one half the curl of the streaming velocity itself [6]. 
However, if the couple force, that is, the force component producing pure rotation, is large, 
the rotation must be treated as an independent variable. The extended description is known 
as Cosserat (or micropolar) continuum mechanics [15, 16], hrst formulated by the Cosserat 
brothers [17, 18] in the late 19th century. Cosserat continuum theory is used in various areas 
such as liquid crystal studies [19] and blood hows [20], and was studied intensively in the 
1950s to 1970s, see Refs. 7, 10, 21-23, 25, and 26. For some reason it is not, however adopted 
by the nanohuidic community. We show that Cosserat theory must also be used for huid 
hows in extremely small conhnements where the molecular structure becomes important. 

(ii) Classical hydrodynamics is based on local constitutive models relating huxes to ther¬ 
modynamic forces. For a shear how the stress at some point depends on the strain-rate 
at that particular point. If the stress depends linearly on the strain rate, this leads to the 
Newtonian law of viscosity [6]. A more general constitutive relation is to let the stress be 
a function in the entire strain rate history and spatial distribution, he., given by a spatial 
and temporal convolution integral of a viscosity kernel and the strain rate [28]. This is the 
approach of generalized linear response theory [3, 29]. The viscosity kernel accounts for the 
characteristic length scale of the spatial correlations [31, 32]; we show below that this must 
be taken into account in order to arrive at the correct huid response on molecular length 
scales. 

Our presentation is based on comparisons of continuum predictions with atomistic molec¬ 
ular dynamics (MD) simulation data. These two descriptions are fundamentally diherent 
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in two ways. First, in MD the system is characterized by discrete particles where the path 
of each individual particle constituting the fluid is traced out through classical mechanics 
[16], he., the particle interactions must be known. The discretization of matter is, of course, 
in strong contrast to the fundamental assumptions of continuum mechanics. Secondly, the 
continuum description applies constitutive relations to form mathematical closed problems. 
No such models are enforced in the standard MD simulations. Any discrepancy between MD 
and the continuum description may therefore be a result of a breakdown of the constitutive 
relation rather than a break down of the continuum theory as such. Our basic conjecture is 
that MD acts as an idealized numerical experiment, and if a given continuum theory agrees 
with the MD data, the theory correctly accounts for the phenomena we study. 

Let us specify, more accurately, what is meant by continuum theory. Basically, one refers 
to deformable fluid volumes characterized by quantities which are continuous at any point r 
over the entire volume and at any time t [6]. This means that these quantities are described 
mathematically by held variables. The basic continuum hypothesis is that one can associate 
a given huid sub-volume (or “huid particle”) with the same characteristic quantities of the 
entire deformable huid volume, no matter how small the sub-volume [2, 6]. Lautrup [7] 
suggests a lower limit of the order of 10^ molecules as stated above, but time averaging 
allows an arbitrarily small huid particle volume as seen in Fig. 1. One held variable is the 
streaming velocity, which is the mass-weighted average velocity of the individual molecules 
in the huid particle around a given point [3]. The huid’s dynamics is governed by balance (or 
conservation) equations. In general the balance equation for some quantity per unit mass, 
0 = 0(r, t), reads in the Eulerian diherential form [7] 

^^ + V- = V-J*. (1) 

where p is the mass density, a,/, a production term, u the streaming velocity, and 3^ the hux 
of 0. Here 0 can be a scalar or vector quantity. The right hand side of Eq. (1) is the sum 
of the body force and the surface force densities, that is, forces per unit volume. When 0 
represent the velocity held, 0 = u, Eq. (1) is the momentum balance equation. The body 
force density can be a gravitational-like force driving the how as in Fig. 1, and the surface 
force density is the pressure tensor, Ju = P [3]. A special case is the mass balance equation 
for which 0 = 1. Since rotation is treated as an independent variable, a balance equation 
on the form of Eq. (1) must be formulated for rotation; this is done in the Supporting 
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Information. (SI) Importantly, in the extended Cosserat description the pressure tensor P 
need not be symmetric [7, 21-23] as in the classical continuum theory. 

Comparison between the continuum description and MD simulation data is carried out 
for molecular fluidic systems at equilibrium, as well as for steady flows in a slit-pore. Here 
we investigate four molecular fluids: a methane fluid, a generic di-atomic (dumbbell) fluid, 
liquid butane, and liquid water. For methane 75 percent of the mass is centered in the 
carbon nucleus and methane is here considered as a simple spherical point-mass molecule, 
as it was done in Fig. 1. Water will, on the other hand, be treated differently using the 
flexible SPC/Fw water model [34] that accounts for the molecular structure and hydrogen 
bonds and thus for the structure of liquid water. The butane model is a coarse grained 
model where the methyl and methylene groups are represented by a united atomic unit, he., 
a spherical point-mass. Details about the butane model can be found in Ref. 35, however, 
here flexible bond are implemented with parameters from the Generalized Amber Force Field 
[36]. The simulations are done using the seplib library [37]. 

Nanofluidic flows are often associated with fluid slippage at the wall boundary [38]. Just 
like the effect of the fluid-fluid interactions on the flow is lumped into a single parameter, 
e.g., viscosity, one effect from the fluid-solid interaction can be modelled into a friction 
coefficient determining the boundary slip. The slippage has a large effect on the flow rate in 
extreme conhnement and is usually quantihed by the slip length Lg. For a Hagen-Poiseuille 
flow in a tube with radius R the relative flow enhancement due to the slip is given as 

[39] 

= 1 + iL./R . (2) 

Lg is typically in the order of a few nanometers. Thus, for a given non-zero slip length the 
flow enhancement increases hyperbolically as the tube radius decreases. The slip is always 
present, but has insignihcant effect on the flow rate for tube radii above microns. Lg is 
normally independent of system size, that is, is not an intrinsic nanofluidic phenomenon 
and is therefore not addressed in this paper. Slippage is here modelled in an ad-hoc fashion 
as it was done in Fig. 1. 

In SI we derive the Cosserat extended continuum theory from the microscopic point 
of view using a microscopic hydrodynamic operator. The derivation, which is based on 
the fundamental dehnition of the macroscopic field variables in terms of the corresponding 
molecular quantities, follows the idea of Irving and Kirkwood [5], and Evans and Morriss 
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[3], see also Ref. 4. The derivation leads to a molecular interpretation of the fluxes entering 
Eq. (1). The hnal dynamical system of equations is sometimes referred to as the extended 
Navier-Stokes (ENS) equations 

= CTJ - Vpeq + + ?7o/3 - rir)V{V ■ u) + {po + ?7^)V^u + 2r]rV x Q, (3a) 

Dr2 

= cTg + 2?7r(V XU — 2S7) + {(y + (q/3 — Cr)V(V ■ Q) + (Co + Cr)V Q , (3b) 

where, D/Df is the material operator, f2 the spin angular velocity held, and Vp is the 
pressure gradient. The transport coefficients rjo, and rjr are the bulk, shear, and rotational 
viscosities, respectively, and Cd, Co and Cr the corresponding spin viscosities. Finally, ctj and 
cTg represent production terms of linear and spin angular momentum, respectively. We refer 
the reader to SI for a derivation and discussion of Eq. (3). 

The theory is only strictly true for isotropic systems, and we study such cases in Sections 
II and III, comparing theoretical predictions with MD simulation data. Flows in extreme 
conhnement are characterized by strong density inhomogeneities and anisotropy. We study 
such hows in Section IV, again comparing theory with MD data. Finally, Section V gives a 
brief summary. 


II. COUPLING: MULTISCALE RELAXATION PHENOMENA IN MOLECULAR 
FLUIDS 

The purpose of this section is to demonstrate the validity of the ENS equations, Eqs. 
(3), by comparing the predictions of diherent thermally induced relaxation phenomena with 
MD simulation data, see also Refs. 42 and 43. We start with this problem instead of 
the situation with conhning walls as the latter is introduces density inhomogeneities and 
molecular alignment at the wall-huid interface. We return to this more complex situation 
in Sect. IV. 

Rather than investigating the quantities directly, one typically studies the associated cor¬ 
relations [44]. Here we will use the approach based on Onsager’s regression hypothesis [45], 
which states that thermal perturbations on average decay according to the deterministic hy¬ 
drodynamic equations of motion. Specihcally, we will compare mechanical spectra obtained 
from MD simulations with predictions from theory. 



A. The stochastic ENS equations 


In equilibrium a fluctuating quantity A can be written as A = Aav + SA, where Aav is the 
average part and SA is the fluctuating part. In equilibrium the average streaming velocity 
and spin angular velocity are both zero so u = (5u and 17 = 50. The fluctuations are mod¬ 
elled using the stochastic forcing approach [46]. Here an uncorrelated zero mean stochastic 
force is added to the constitutive relations, see SI. For example, for the antisymmetric pres- 

ad ad 

sure the constitutive relation with stochastic forcing reads P= —r]r{'V x 5u — 250) -|- 5 P, 

ad 

where 5 P is the fluctuating part of the flux. 

To a first order approximation in the fluctuation we have on the left-hand side of Eq. (3) 


, ^ ,D5u (95u 

{Pav + ^ Pav-^ 


( 4 ) 


and 


{Pav A 5p) (^lav A 5/) 


D50 

i5r 


Pav I a 


950 


( 5 ) 


In Fourier space the stochastic ENS equations read to first order in the fluctuations for wave 
vector k 


Pa^~^ = -*kpeg - iVv A ?7o/ 3 - 77r)k(k ■ 5u) - (? 7 o -h rir)k‘^6u + 2ir]rk x 50 -h fk ■ 5P 

(6a) 

95 O - - - - -- r, -- -- Cld 

=2i;,,(!k X (5u - 26(1) - ((„ + (o/3 - (r)k(k ■ Sft) - ((o + Qk S(l + ik ■ iQ + P , 

(6b) 


due to the properties of the divergence operator. See SI Eq. (2) for the definition of the 
Fourier transform. It is here convenient to introduce the following coefficients 


Vt = VoA T]r, 0 = Co + Cr, and Cz = C« + 4Co/3 , (7) 

where subscripts t indicates “transverse” and / indicates “longitudinal”. It has been shown 
[10, 42] that Ct ^ Cb and we write both coefficients (. We may then define the susceptibility 

Xik"^) = 4r/r + , (8) 

where = k^. We will also drop the subscript av from here on. 
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We take k = (0, 0) and write out the x-component of the velocity and 2 ;-coniponent 

of the angular velocity 


P 

pi 


dSu 

d6n 


X 


z 


—rjtk‘^5ux + 2irjrky5VLz + ikySPyx 
-x{k‘^)SQz - 2irjrky5ux + ikySQy^ + 2S . 


(9a) 

(9b) 


These two components are both transverse components to the wave vector and are coupled. 
We also investigate the longitudinal angular velocity component SQy which is given through 

d5kl,, n -- -- Cld 

pi = ~xik )6kly + ikySQyy + 25 Py. (10) 


Note, that this longitudinal component is unaffected by the coupling between the linear and 
spin angular momenta. 

We dehne the following three correlation functions 


CL(k,t) = (i5Mi(k,«)i5Mi(-k,0))/r 

(11a) 

Cs)-„(k, 4 ) = («j,(k.«)fa,(-k.0))/r 

(11b) 

ci;„(k,() = {«i,(-k,()S„(-k,0))/r, 

(11c) 


which we denote the transverse velocity auto-correlation function (TVACF), the transverse 
cross-correlation function (TCCF), and longitudinal angular velocity auto-correlation func¬ 
tion (LAVACF), respectively. By assumption the fluctuating fluxes are uncorrelated with 
the velocity and angular velocity, e.g., (SPy^Su^) = 0. Thus, multiplying Eqs. (9a) and Eq. 
(9b) with (5Ma;(—k, 0) and ensemble averaging we arrive at the differential equation system 
for the TVACF and the TCCF 

fif-L 

P—^ = -Vtk'^CL + ‘liPrkyC^^ (12a) 

fif-L 

- ‘liVrkyC^^ . (12b) 

Similarly, multiplying Eq. (10) with (5r2j^(—k, 0) one has for the LAVACF 

pl^ = -x(k^)C'an (13) 

upon ensemble averaging. Now, Eqs. (12) and (13) can be solved yielding to second order 
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in wave vector 


cLiKt) = 

{Ik^{ 

(14a) 

cLiKt) = 

i2prk 

(14b) 

Apr + {I{Vr-Vo) + Ok^ 


QksT ^ 

Apl 

(14c) 


where the characteristic frequencies are 

+ IVrk"^ j Vok^ 

ojq = ——, uJi = --, and u :2 = - 

pi pi p 

The pre-factors in Eqs. (14c) and (14a) are calculated from the hrst order approximation 
in the fluctuations J ~ p(5u as above and therefore 



(5u(k, f) ^ ^ CjC (16) 

k i 

from the dehnition of the linear momentum density, see SI Eq. (14) in SI, and for one 
component systems with molecular mass m. Likewise, to a hrst order approximation in 
density and moment of inertia pS ~ p/Li, and from SI Eq. (29) 

(17) 

i 

as / = 2Ip/3 [11]. Applying the equipartition theorem one arrives at the pre-factors. Equa¬ 
tions (16) and (17) also provide a hrst order method to calculate the correlation functions 
in the MD simulations; this method is used here. 

It is informative to work in the frequency domain, i.e., to predict the peak frequencies in 
the corresponding spectra. Applying the Fourier-Laplace transform dehned by 


we get 


poo 

£[/](k.w)=/ /(k,«)e-"‘d« 

Jo 


Cuu%^) 

C^on(k,ca) 


_kBT f Ik^ _4 + /A;2\ 

4p ycai + iuj U2 + iui J 

i2r]rk f ^ ^ ^ 

x(A;2) j(^r]r - r]o)kJ \a;i + iu U2 + iuj J ' 

Qkj^T 1 
Apl uq + ioj 


(18) 


(19a) 

(19b) 

(19c) 
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From Eq. (19b) we can make a very important conclusion, namely, 


Cn«(k,ii') ^ 0 for k ^ 0 , 


(20) 


This means that the coupling can be ignored on long length scales. This is also expected 
as the classical Navier-Stokes theory holds for macroscopic systems and no coupling effect 
is observed. The relaxation of spin is still governed by the rotational viscosity, but this 
relaxation does not affect the relaxation of linear momentum controlled by the usual viscous 
dissipation processes. If we define Uc as 


Wc = hmcuo = —, 
k^>0 pi 

the LAVACF and TVACF are, in the limit of zero wave vector. 


( 21 ) 


= C^kw) = ——(k-»0), (22) 

Apl UJc + toj p U2 + toj 

Furthermore, for the fluids studied here the effect of the coupling on the TVACF is not 
large, that is. 


Ik^ 


oji + ioj 


< 


4: + le 


0J2 + ioJ 


(23) 


even for wave vectors in the sub-molecular diameter range, and the limit in Eq. (22) need 
not to be taken as a strict limit. It is also worth noticing that the rotational viscosity pr is a 
linear function of the moment of inertia I for sufficiently large I [17, 43], so Wc is independent 
of I here. 


B. Comparison with molecular dynamics 

We hrst compare the predictions from the continuum ENS theory with MD simulation 
data for the simple di-atomic molecule (the dumbbell model) in the super-critical fluid 
regime. The transport coefficients, po, and pr are listed in Table I in SI. 

Figure 2 shows MD data (symbols connected with lines) for imaginary parts of the spectra 
of the TVACF and LAVACF; normalization is carried out for clarity. The prediction from 
the continuum theory is plotted as full blue lines. It is observed that for small wave vectors 
k = 27r/L the continuum prediction is in excellent agreement with the MD data, but it 
fails for larger wave vectors k = 207r/L. We emphasize that no htting is performed, and 
all relevant parameters are taken from SI Table I found from independent simulations and 
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(a) 


(b) 



FIG. 2. (Color online) Spectra for the dumbbell model. Blue lines are predictions from the con¬ 
tinuum theory (CT) using the coefficients given in Table I, SI. The arrows indicate peak frequency 
behavior with increasing wave vector, (a) The imaginary part of the spectrum for the longitudinal 
angular velocity auto-correlation function (LAVACF). Normalization is carried out for clarity in 
the comparison. The dashed line indicates ojc given by Eq. (21) (b) As in (a), but for the transverse 
velocity auto-correlation function (TVACF). The dashed line indicates U 2 given by Eq. (15) (c) As 
in (a), but for the transverse cross correlation function (TCCF). The predictions from the theory, 
Eq. (14b), for small wave vectors k < 27r/L are shown. Typical orders of magnitude for the MD 
units are cr = 1 A, m = 10“^® kg and e/ks = 10^ K. L = 13.17a. 
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methods. Using typical values for the MD units a, e and m the results show that the 
continuum theory predicts the mechanical spectrum for wave lengths in order of 2-3 nm and 
above, and time scales in order of of 1-10 ps and above. 

For the TVACF, Fig. 2 (b), the result can be understood from the fluid stress relaxation 
at zero wave vector as suggested by Bocquet and Charlaix [2]. From the last equation in 
Eq. (15) we can define a wave vector dependent relaxation time T 2 = 2np/. This 
relaxation time must be larger than the characteristic relaxation time r* at zero wave vector 
for the predictions to hold < T 2 , he., for the viscosity to be wave vector independent. 
This means that 


VoTsk'^ 


<1 or k < i I . (24) 

27rp V VoTs 

We will denote this the Bocquet-Charlaix criterion. Estimates for the relaxation time is 
given through the shear pressure (or equivalently stress) auto-correlator 


OS OS 

G(t) = (t) (0)) , 


(25) 


where P^y is the {x, y) component of the symmetric part of the pressure tensor. For the 
dumbbell model G{t) is fully decayed at Tg ~ 3a/\/m/e which gives k < 1.2a~^. This is in 
perfect agreement with the results depicted in Fig. 2 (b). Alternatively, the relaxation time 
can be given through the Maxwell relaxation time tm = Po/G°^ or the viscous relaxation time 
[49] Ty = 4/i o/277o, where G°° is the inhnite shear modulus and 4/i q is the hrst normal stress 
coefficient. For the diatomic model studied here tm ~ Ty = 0.05(TA/m/e giving k < 9Aa~^ 
which is not what is observed. Therefore, the characteristic decay time that should be used 
for the Bocquet-Charlaix criterion is the time for the autocorrelation function G{t) to fully 
decay. 

In the small wave vector regime the relaxation of spin angular momentum is dominated 
by the coupling mechanism between linear and angular momenta as the peak is located at 
c^o ~ = 4:Vr/{pI)- The relaxation of linear momentum. Fig. 2 (b), is on the other hand 
due to usual viscous mechanisms seen by the peak frequency U 2 = rjok'^/p. For large k 
the continuum theory overestimates, by an order of magnitude, the peak frequency for the 
LAVACF and TVACF due to over-estimation of the effect of the spin diffusion. 

Figure 2 (c) depicts the TCCF for the dumbbell model. Again, the theory performs 
surprisingly well for sufficiently small wave vectors {k < 2n/L), but fails for larger. It is 
worth noticing that the amplitude of the TCCF is a non-monotonic function with respect 
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FIG. 3. (Color online.) Liquid butane. Molecular dynamics results for (a) the LAVACF spectrum, 
and for (b) the TVACF spectrum. The dashed lines indicate Uc, Eq. (21), and UI 2 , Eq. (15). The 
arrow indicates peak frequency behavior with increasing wave vector. L = 32.95 A. 

to wave vector, having a maximum around k = Att/L. This behavior is also captured by 
the ENS theory. To illustrate that the amplitude is a decreasing function of wave vector the 
TCCF for k = 71/L and k = 7r/(2L) is plotted as predicted by the theory. Recall, in the 
limit k —)■ 0 the coupling vanishes. 

Next we apply the theory to liquid butane. As discussed above the butane model is not 
uni-axial or rigid, however, from the principal moment of inertia we argued that the theory 
should be a good approximation. The result is shown in Fig. 3. For the LAVACF one 
observes a peak frequency at around 9 THz and the relaxation process is extremely fast. 
This fast mode is not precisely captured by the theory with Wc = 6.6 THz. Interestingly, 
the peak frequency is almost independent of wave vector for the range studied here. This 
indicates that for these fast modes the diffusion of spin is less important for the relaxation 
processes. For the slower relaxation of linear momentum, we see that the peak frequency 
is predicted very well by the theory, the MD result is a;=1.0 THz and the predicted one 
is around U 2 = 0.92 THz for the lowest wave vector. Again, for larger wave vectors the 
prediction fails as expected. 


III. NON-LOCAL RESPONSE 

The classical linear constitutive relations are local in the sense that the flux only depends 
on the local and instantaneous thermodynamic force. This is in general not the case, rather 
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the response depends on the entire force distribution in the system, as well as its history. 
One can model this phenomenologically by introducing frequency and wave vector dependent 
transport coefficients [3]. This renders the continuum description valid on arbitrary small 
length and time scales. The generalized transport coefficients are referred to as kernels. In 
the homogeneous isotropic case, assuming space and time invariance, the linear non-local 
constitutive relation for the symmetric part of the pressure tensor reads [3, 50] 

os /■* 

P(r,t) = — / / ? 7 (r — r', t — t') 7 (r', t)dr'dt' . (26) 


' —CX) O —OO 


7 (r, t) =Vu (r, t) is the symmetric part of the velocity gradient, he., the strain rate. Fourier 
transforming with respect to space and Fourier-Laplace transforming with respect to time 
yields 

05 ' 

P(k, u) = -r/(k, a;) 7 (k, u) (27) 


from the convolution theorem. 

The shear viscosity kernel I/(k, ca) can be found from the TVACF as it is now shown. 
Here we focus on molecules with small moment of inertia and small wave vector regime, he., 
small Ik^. In this limit Eq. (19a) can be rearranged giving 


h(k,w) 


kpT -iupC^S\^,u) 


(small J/c^) . 


(28) 


In particular, we have at zero frequency 


^ (small JA;^). (29) 

This approximation holds even for large values of as discussed above, Eq. (23). For 
molecules that can be regarded as point masses, say methane, the moment of inertia is zero 
and Eq. (28) is exact. In Fig. 4 (a) the viscosity kernel at zero frequency is plotted for the 
methane, dumbbell, butane, and water systems using Eq. (29). One immediately notices 
that for k ~ l(j~^ the wave vector dependent viscosity approaches the zero wave vector limit. 
This is in good agreement with the Bocquet-Charlaix criterion, Eq. (24). Interestingly, this 
is independent of the specihc fluid studied here and the local constitutive relations can be 
applied on length scales down to approximately 27r//c ~ 2-2.5 nm. 

Is this a general result that applies to all fluidic systems? The answer is no! In Fig. 
4 (b) the zero-frequency viscosity kernel is plotted for the asymmetric dumbbell model for 
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FIG. 4. (Color online) (a) Viscosity kernels for the dumbbell model, butane, water, and methane. 
For butane, water and methane cr = 3.9233 A, 3.166 A and 3.80 A, respectively, (b) Viscosity 
kernel for the asymmetric dumbbell model at different temperatures. The Newtonian viscosity 
values are tjq = for T = 2.QelkB and rjo = ASi^fraija^ for T = For (a) and 

(b) the dashed lines are best fit to the empirical form rj{k) = r/o/(l + ak^) [51] and is included to 
guide the eye. 


different temperatures. The asymmetry arises due to the mass and Lennard-Jones parameter 
differences between the two constituent atoms. The asymmetric dumbbell model allows 
one to probe the dynamics in the highly viscous regime without crystallization occurring 
[52]. The result shows that for relatively high temperatures the kernel has the same wave 
vector dependency, but approaching the viscous regime (lower temperature) the kernel only 
reaches the Newtonian viscosity at longer length scales. This indicates that the dynamical 
processes behind the viscous response take place on longer length scales in accordance with 
the cooperative motion in super-cooled liquids [53]. The non-local viscous response has also 
been studied for highly viscous two component Lennard-Jones system and polymer melts, 
see Refs. 31 and 32. 

The failure of the local constitutive relation, that is, of Newton’s law of viscosity, is very 
clearly illustrated by Todd et al. [54] for a point-mass Weeks-Chandler-Andersen (WCA) 
system [55]. In real space the non-local description amounts to a convolution of the viscosity 
kernel and the strain rate distribution, Eq. (26). In the homogeneous situation where the 
fluid undergoes a steady shear in the a:-direction with varying amplitude in the ^-direction 
we have one non-zero shear component in the pressure tensor, namely, the Pxz component. 
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In this steady situation Eq. (26) reduces to 


Pxz{z) = - rj{z- z')j(z') dz' , (30) 

J —OO 

where 7 ( 2 :) = dux(z) / dz. If the shear is induced by an external force held Fe{z) = Fq cos{kz), 
the huid how is Ux{z) = cos{kz), where is the excited Fourier mode of the velocity held. 
We assume that this is the only mode excited, he., the force amplitude must be sufficiently 
low [56]. Also, this ensures a linear response as well as constant temperature and density. 
The strain rate is then 

7 ( 2 ;) = —ku^^sinikz) . (31) 

For simplicity we shall assume that the kernel is given by a Gaussian function 

V{z) = (32) 

such that 1/y/a gives a characteristic decay length. The kernel must fulhll [54] (i) 
f^7](z) dz = r/o, and (ii) ?](z) is an even function. Substituting Eqs. (31) and (32) 
into Eq. (30) we have upon integration 

Pxz{z) = r]Qku^e~^^sm.{kz) . (33) 

If rj{z) = rjoSi^z) the model is local, corresponding to Newton’s law of viscosity, that is, for 
the local model 

PL{^) = -Voi = mkul sin{kz) . (34) 

The system can be simulated using the Sinusoidal Transverse Force (STF) method [57], 
and it is possible to evaluate for diherent external force helds and wave vectors. The 
two diherent predictions can be compared to the actual shear pressure P^, which is found 
directly from the momentum balance equation that for the steady how reads 

= pP. . (35) 

Integrating we obtain 

Pxziz) = ^ sm{kz) . (36) 

The comparison is made in Fig. 5. Clearly, the local prediction fails for the larger wave 
vector. Fig. 5 (b), whereas the non-local prediction agrees with the actual shear pressure. 
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FIG. 5. (Color online) Pressure profiles for the Weeks-Chandler-Andersen point-mass system under 
periodic shearing force at state point (p, T) = (0.685cr“^, 0.765e/A:B). a = 4.81cr“^ is found from 
best fit to data given in Ref. [56]. Input data for the theory is taken from Todd et al. [54]. 
(a) Small wave vector: k = 0.357cj“^,= 0.887^/e/ m,Fo = 0.15e/{am) (b) Large wave vector: 
k = 3.57a~^ ,u^ = 0.027s/ejm, Fq = 0.225e/(am) 

From the non-local model we conclude that spatial correlations result in a reduced shear 
pressure. 

From Eqs. (33) and (34) we can quantitatively evaluate the effect of spatial correlations 
on the stress. Specihcally, we have the relative difference given by 

AP^f = 1-^ = 1- . (37) 

^ XZ 

For the WCA system studied here ~ 0.0a^Jm/e and the Bocquet-Chairlaix criterion gives 
k < 2.8cr“^, this corresponds to an error in the stress below 32 % according to Eq. (37). 

The Gaussian function does not perfectly fit to the kernel data. Nevertheless, this simple 
functional form captures the non-local response well, due to the smoothing of the convolu¬ 
tion. Other more complicated forms have been suggested, see Refs. 31, 51, and 58. 

Todd and Hansen [59] showed that the non-local response is only relevant for flows where 
the strain-rate is non-linear with respect to position. Couette and Poiseuille flows are then 
not affected by non-locality. To illustrate this consider any functional form for the kernel 
which fulfills the criteria given above: its integral gives the zero wave vector viscosity and it 
is an even function with respect to 2 ;. First, making the change of variables u = z — z' Eq. 
(30) reads 

poo 

Pxz{z) = - r]{u)'y{z - u)du. (38) 
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Assuming a strain rate on the form j(z) = az, we have 


PrAz) = —az 


T]{u)du + a / r]{u)udu =—ar]oz , 


(39) 


as the integrand in the second integral is an odd function. This result is the same as the local 
prediction. In general, if a Taylor expansion of the strain rate 7 = oq + aiz + ... + a„ 2 :” + ... 
exists, using the properties of odd and even functions one can verify that the non-zero 
non-local effects of the strain rate can be determined by the even moments of the kernel [59] 


Mn = / rj{z)z^ dz (n > 0 and even) 


(40) 


In the case of Couette and Poiseuille flows the Taylor expansion terminates at zero’th and 
hrst order, respectively, and there are no non-local effects. 

The spin and rotational viscosity kernels can be found by simply rearranging Eq. (19c) 
giving the generalized susceptibility 


X(k,cu) 


C|ln(k,i..) 


(41) 


Our group recently [42, 43] conjectured that the rotational viscosity governs the fast wave 
vector independent relaxation processes as indicated in Figs. 2 (a) and 3 (a). This transport 
coefficient is therefore only frequency dependent. We then have 


x(k,w) = 4))r(w) + ((k.w)*:^ 


(42) 


and therefore 


rjriu) 


--hmY(k,a;) and C(k,a; 

4 k^O ^ ^ ^ 


X(k,a;) -driAuj) 


(43) 


We called this the generalized extended Navier-Stokes (GENS) theory. From MD simulations 
one can calculate the LAVAGE (as shown above) and from there hnd the kernels. For dense 
fluids ( is characterized by a sharp peak around zero wave vector [42] since the diffusive 
contribution to the relaxation of the LAVAGE is very small for k > see Figs. 2 (a) 

and 3 (a). The spin viscosity kernel has the same properties as the shear viscosity kernel 
and for this reason we do not expect any non-local effects for flows where the gradient of 
the angular velocity is constant or linear. 


20 



IV. NANOFLOWS 


A. The Poiseuille flow 

We first study a Poiseuille flow, the geometry is shown in Fig. 1. In experiments this 
flow can be achieved by application of a constant pressure gradient. Generating a pressure 
difference in simulations with, for example, a piston and using molecular reservoirs can cause 
density variations in the direction of the flow and other inlet/outlet effects. We therefore 
use a constant force held acting on each point mass in the fluid to drive the flow. The wall 
particles are arranged on a simple cubic lattice and are allowed to vibrate around their initial 
lattice site using a simple restoring spring force. The viscous heating generated in the fluid 
is removed by thermostating the wall particles. This method resembles the real physical 
experiment and is therefore often referred to as direct non-equilibrium molecular dynamics. 
The interested reader is referred to Ref. [11] for further details. To get a satisfactory signal- 
to-noise ratio in the MD simulations unrealistically large external forces are typically applied 
to drive the system, and the resulting flow rates are very large, typically, in the order of 
10-100 m s“^. Despite these large flow rates the Reynolds number is usually less than unity 
due to the extremely small characteristic length scales involved. Finally, it is very important 
to ensure that the simulations are carried out in the linear regime which is discussed below. 


1. Continuum predictions 

In the linear regime of low Reynolds number and for the geometry shown in Fig. 1 the 
ENS equations form a two-point boundary value problem in the steady state 

^ d?Ux dD^ /. . , 

pFe + ht-rw ~ ^F-r- = 0 (44a) 

dz^ dz 

2,,(!^_2fi,)+C^ = 0. (44b) 

where —h < z < h. Recall that pt = Vo + Vr and ( = Co + Cr- Introducing z' = z/h, —l<z'< 
1 and applying no-slip boundary conditions, Ux{—1 ) = u^iX) = 0 and 1) = Dy(l) = 0, 
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Eringen [23] solved this yielding 


uJz') = Uci I - z + 


,2 2rirCoth.{Kh) f cosh{Khz') 


rjtKh \ cosh.{Kh) 


O ( — f _ , 

^ ^ h \ sinh(i^/i) 

with the following dehnitions of Uc and K 


Uc = h‘^pFe/{2riQ) and K = {^VrVo/iCVt)) ^ 


(45a) 

(45b) 


(46) 


The application of the no-slip boundary condition is not justihed. A correct treatment ap¬ 
plies the Neumann boundary condition for both the velocity and angular velocity helds, 
however, this is not straightforward in that the two are likely coupled and therefore depen¬ 
dent on each other. While the boundary condition for the velocity held has been studied in 
great detail, see e.g. Refs. [60-62], very little is known about the spin boundary condition. 
Just recently De Luca et al. [63] showed that the spin held does possess slippage and Badur 
et al. [64] used spin slip to account for how enhancement. As mentioned in the introduction, 
we will treat the problem in an ad hoc fashion and simply set the angular velocity slip in 
accordance with the MD data. 

If one ignores the coupling, rjr = 0, the solution for the streaming velocity, Eq. (45a), 
reduces to the classical Poiseuille how solution 


u. 


,(z') = !lc (l - . 


(47) 


In this classical situation the angular velocity is found from the vorticity, VLy = ]^dux/dz, 
that is, 

(48) 


Cl / r\ / hpFf, ! 

^ h 2rio 


in agreement with Eq. (45b) for pr 0. As the classical treatment does not allow for 
specihcation of spin boundary condition, Eqs. (45b) and (48) diher by a magnitude of 
hpFe/{2po) at the walls. 

From Eq. (45a) one can see that the maximum velocity, located at z' = 0, is lowered 
as a result of the coupling since from the last term we have 1/cosh(iLh) — 1 < 0 and thus 
Mj,(0) < Uc- Another way to quantify this ehect is to evaluate the volumetric how rate Q [3] 


/ w ph ph 

/ u^{z)AzAy ^ 2w / u^iz)dz , 


(49) 


' —w J —h 


f-h 
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FIG. 6. (Color online) Relative volumetric flow rate reduction for the dumbbell fluid, liquid butane, 
and liquid water, a = 3.92 A and 3.17 A for butane and water, respectively. 


where w is the characteristic half length in the ^/-direction. This gives the relative volnmetric 
flow rate rednction 





3r]r(tanh.{Kh) — Kh) 
Tjt tanh.{Kh){Khy 


(50) 


Equation (50) is plotted in Fig. 6 for the dumbbell fluid, liquid butane, and water. The 
relevant coefficients can be found in SI Table I. It can be seen that for water flowing in a 
channel with a width of 9 nm the flow rate is reduced by about 10 % due to the coupling. 
As the channel width increases the flow rate approaches that of the classical predictions and 
the effect of the coupling can be ignored. 

From Eq. (50) the relative flow rate reduction increases as the product Kh decreases. 
From this observation, one can dehne a characteristic fluid length scale Ic [65] below which 
the effect of the coupling becomes signihcant. To this end we write the parameter K as 


K = 


L 


with Ic ^ t/cTT-- 


( 61 ) 


From SI Table I it is seen that po > hr and K k, 2/Ic- Thus a signihcant how-rate reduction 
occurs for huids with large critical length scale Ic- For water Ic = 3.5 nm and for butane 
Ic = 0.5 nm in agreement with the relative large how rate reduction observed for water. 


2. Comparison with molecular dynamics simulations 

First we compare MD data with continuum predictions for the dumbbell system. Be¬ 
fore the comparison, however, the linear Newtonian response regime should be identihed, 
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at least, for the bulk fluid region. To this end one can apply the synthetic SLLOD al¬ 
gorithm developed by Evans and Morris [ 66 ]. Basically the method imposes a constant 
strain rate (linear velocity prohle) on the system while ensuring a homogeneous density 
and iso-kinetic temperature. To achieve this the equations of motion are reformulated ac¬ 
cording to the Gaussian principle of least constraint, see also Refs. 3 and 4. Performing 
a series of SLLOD simulations it is found that the Newtonian regime occurs in the range 
0 < 7 < 0.05((TA/m/e)“^ for the dumbbell fluid . The upper limit for the external force 
field can then be approximated by rewriting Eq. (48) to 7 ~ 2VLy = —2uclhz' giving 
Ff. < 2?7o7m/ph, with 7^ = e)~^. Note, in the wall-fluid region the velocity may 

feature rapid changes and here the linearity is not guaranteed. 

Based on the SLLOD approach Delhommelle [67] developed a synthetic method to cal¬ 
culate the rotational viscosity rjr as a fnnction of spin angnlar velocity. See also Edberg et 
al. [ 68 ]. To onr knowledge no synthetic or controlled method exists to stndy the spin viscos¬ 
ity dependency of the gradient of the spin, or the rotational viscosity dependency on strain 
rate. We will therefore here assnme that the linear regime is identical to the Newtonian 
regime, i.e., in the regime where the viscosity is independent of the strain rate. 

The classical description predicts that the Poisenille flow is local flow according to Sec. 
III. Also, from Fig. 6 we expect the flow-rate redaction dne to the conpling to be very low for 
the dnmbbell model. Thus, based on the theory we can expect the classical description to be 
a good approximation for this system. The time-averaged velocity and spin angular velocity 
prohles are shown in Fig. 7 (a) for the dnmbbell model where the pore width is approximately 
14.8 atomic diameters. The profiles are sampled after the system has reached the steady 
state. The temperature prohle (not shown) is constant and the temperature is T = d.Oe/Zcs 
throughout the channel. The predictions from the classical Navier-Stokes theory, Eqs. (47) 
and (48), are also plotted using the shear viscosity from SI Table 1. Velocity slippage at 
the wall-hnid interface is allowed, Ux(z) = Ucil — {z/hY) -|- is then the only htting 

parameter in the comparison. The agreement between MD data and classical continnum 
predictions is excellent, except at the wall-hnid interface. This is highlighted by the shear 

OS 

pressnres plotted in Fig. 7 (b). According to the classical theory the shear pressnre Pxz 
is linear, however, at the wall-hnid interface this is not the case. Also, the classical theory 

ad 

assnmes a zero anti-symmetric part of the shear stress Py. This is clearly not fnlhlled near 
the wall. 
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FIG. 7. (Color online) (a) Velocity and angular velocity profiles for the dumbbell model undergoing 
a Poiseuille flow. Symbols represent MD data, and lines the classical predictions where slippage is 
included, (b) The corresponding shear pressures. Pxz is calculated by integration of the momentum 
balance equation, using the density profile given in (c). The antisymmetric pressure is calculated 
from the constitutive relation SI Eq. (40) using MD data as input, (c) Density and order parameter 
profiles. The molecules in the lower left corner illustrate (exaggerated) the molecular ordering near 
the wall. 


To understand the disagreement at the boundary, we analyze the fluid ordering. It is 
well known that the wall induces a density variation in the fluid [69]. The density prohle is 
shown in Fig. 7 (c) (black dots). It is seen that the density varies in a region approximately 
one atomic diameter away from the wall. The transport properties are functions of density, 
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and one should expect a variation in the viscosities here. Furthermore, one can evaluate the 
molecular alignment ordering through the parameter [70] 

p=^cos\e)-^, (52) 

where 6 is the angle between the molecular bond and {x, |/)-plane. For perfect parallel 
alignment p = —1/2, that is, the molecules closest to the wall are, on average, aligned with 
the wall. For distances around one atomic diameter the molecules are slightly normal to 
the wall as p > 0. The extremes are illustrated with the two molecules in the lower left 
corner in Fig. 7(c). For zero order parameter, the molecules have random orientations 
which is the case in the interior of the channel. This means that the system possesses a 
degree of anisotropy in the wall-fluid region. To fully account for the density variation and 
ordering one should therefore describe the transport properties through a position-dependent 
tensorial shear viscosity. 

Figure 8 (a) shows the velocity and density prohles for a butane flow where the pore width 
is just 6 nm. For such extreme conhnements the fluid layering stretches over the entire pore. 
The order parameter prohle. Fig. 8 (b), shows that the molecular orientation is strongly 
anisotropic. Finally, the mean square molecular end-to-end distance also varies, showing 
that the butane molecule on average is elongated at the fluid-region by around 2 per cent. 
Such a complex system is not modeled appropriately by the classical or extended theories 
presented here. This is not an indication of a breakdown of the continuum picture, but an 
incomplete modeling. Worth noting is that the fluid ordering and layering is constant over a 
large range of external forces including zero force, see Fig. 8, and is thus not flow induced. 

As pointed out by Bitsanis et al. [71] the velocity prohle features surprisingly small 
modulations considering the density prohle: one should expect the transport properties to 
vary signihcantly across the channel having large ehects on the how prohle. The authors 
suggested the local average density model (LADM) wherein the transport properties at a 
point is a function of the average density around that point. In the current geometry, 
where the density is constant in the plane parallel to the wall, the local average is 

'pi.z) = -r / p{z')dz' , (53) 

A A-Ia 

where A dehnes the region of averaging. The agreement between the LADM and simulation 
data can be very good, especially if one introduces a non-uniform weighting function [72]. 
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FIG. 8. (Color online.) (a) Density and velocity profiles for butane in a slit-pore of width 6 nm. 
The prediction from the theory breaks down and is not shown, (b) Corresponding order parameter 


and square end-to-end distance profiles. 

However, the LADM cannot predict the shear pressure response in Fig. 5 as the density 
is constant. Also, the LADM model is not capable of predicting the strain rate reversal 
observed by Travis et al. [11, 73] 

To account for the observed velocity prohle one can write the position-dependent (in¬ 
homogeneous) non-local constitutive model as 





( 54 ) 


The position dependency likely comes from the varying density at the wall-fluid region. The 
application of this relation is not straightforward [74, 75] as it is unclear how the convolution 
should be performed at the wall where the support of the kernel goes beyond the boundary 
and is unknown [74, 75]. Recently, Dalton et a/. [76] used a sinusoidal longitudinal force 
(SLF), also introduced by Hoang and Galliero [72], to control the density variation in a 
periodic system. Due to the periodicity, the boundary problem can be eliminated. The 
density prohle can be controlled to such an extent that it resembles that seen in conhned 
systems. The huid can then be driven by an STF. The authors showed that the non-local 
response is capable of predicting the strain rate reversal observed by Travis et al. [11, 73], 
as well as the relative small modulation on the velocity prohle. A rigorous and general 
implementation of Eq. (54) into the balance equations for conhned systems is still lacking. 

For the dumbbell and butane models the coupling between the linear and spin angular 
velocities has little ehect on the how. From Fig. 6, however, the ehect is signihcant for 
water how in channel with widths below 5 nm. In Fig. 9 (a) MD data for the velocity 
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FIG. 9. (Color online.) (a) Velocity profiles of water undergoing a Poiseuille flow, (b) Relative 


profile curvature difference at z = 0. From Ref. 9 with modifications. 

profile for water is plotted where the channel width is approximately 10 water diameters. 
Also, shown are the predictions from the classical NS and ENS theories. Density and order 


parameter prohles (not plotted here) show little density variation and molecular alignment, 
except within 3-4 A of the wall. The slip velocities are estimated by htting a second-order 


polynomial (dashed lines) to the velocity prohle, excluding the wall-fluid region where the 
fluid is slightly anisotropic and inhomogeneous; this then amounts to the apparent slip length 
[77] and is the only htting parameter used in the comparison. It is seen that the classical 
prediction fails, while the ENS theory captures the how prohle Note, ht of the prohle data 


to the classical description will result in a wrong viscosity no matter how many data points 


in the wall-huid region are included. Any shift of the prohle will not change this either. An 
extra source of dissipation must be present. Furthermore, note that a complete description 
involves a position dependent non-local anisotropic modelling of the wall-huid region. 

To remove any ehect in the boundary region one can evaluate the curvature in the channel 
midpoint, z = 0. The predictions are simply found from the second-order derivative of Eqs. 
(45a) and (47). The relative diherence is 


T]r coth.{Kh)Kh 



(55) 


(hr + ho) cosh(A'h) 


is plotted in Fig. 9 (b) together with the results from the MD simulations. Within 
statistical uncertainty the ENS theory and MD simulation results agree. As the channel 
width increases the relative curvature diherence vanishes and the classical description is 
re-captured. 

The particular model applied is parameterized with respect to the liquid state and the 
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wall is a Lennard-Jones cubic lattice, see Ref. 9. The fluid structure near the wall and its 
effect on the dynamics will be affected by the different models, choice of model parameters, 
and wall details. However, it is not the aim here to critically review the fluid structure at 
near the wall, but to see the effect of the coupling. 


B. Inserting torqne 


Perhaps the most clear illustration of the translational-rotational coupling is seen by 
introducing an external torque into the system while having a zero production term for the 
linear momentum. In general, if the resulting torque density pT^. is sufficiently small then 
for the geometry in Fig. 1 we have 


pFe -f 2rir 



d^Ux ^ dHj^ 

ht , 2 ^Vr .. =0 

dz^ dz 

(56a) 


(56b) 


Integrating Eq. (56a) we get du^/dz in terms of VLz which is substituted into Eq. (56b) 
resulting in a second order inhomogeneous differential equation for From this and Eq. 
(56a) and by application of Dirichlet no-slip boundary conditions one has 


Ux{z') 

VLy{z') 


dr]r 


Cl sinh.{KhC) 


VtK 

= 2Ci cosh.{Khz') 


P^Vt 


2rir{p^Vt + 2rjrCo)h _ Coh 
2tij.Cq 


CVtK^ 


(57a) 

(57b) 


where Cq and Ci are integration constants. One can show that Ci goes rapidly to zero as h 
increases. In this limit the spin angular velocity is 

o r b = _ dpYeTllrjtK 

and velocity prohle is linear with a slope given by the last term in Eq. (57a). Figure 10 
depicts the two prohles for the butane liquid using Fe = 413 m^s“^. From this one sees 
that the external torque produces a signihcant local flow; the average flow is zero due to the 
system symmetry. 

In 2009 Bonthuis et al. [78] showed that the coupling between the linear momentum and 
spin could be exploited in order to pump water through carbon nanotubes by application of 
a rotating held. The was theory based on the ENS equations at it was noted that in order to 
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FIG. 10. Flow of butane as a result applying an external torque with h =. The insert shows the 
corresponding angular velocity normalized with respect to the value in the limit of large h, Eq. 
(58), here denoted flm- 

obtain a non-zero mean flow asymmetric boundary conditions must be employed which can 
be achieved by conhning the fluid between two walls with different hydrophobicity in the 
case of water pumping. Recently De Luca et al. [63] performed extensive MD simulations 
of the mechanism under experimentally feasible conditions, indicating that the mechanism 
is functional. This could prove to be a way to overcome the large hydraulic resistance 
characterizing nanofluidic flows. Felderhof showed in 2011 that the coupling can also be 
utilized to perform plane-wave pumping [79] and even propel microrobots [80]. 


V. SUMMARY 

We have derived the relevant dynamical equations for isotropic nanofluidic flows. The 
formulation is based on the basic dehnition of a macroscopic field variable from the corre¬ 
sponding microscopic or molecular variable, and it includes the underlying molecular struc¬ 
ture. Two intrinsic nanofluidic phenomena were discussed, namely, (i) the coupling between 
the spin angular momentum and (ii) the linear momentum and the non-local fluid response. 
The important points are the following. 

(1) The effect of the coupling between the linear and spin angular velocities can be 
estimated through the characteristic length scale Ic, Eq. (51). For large Ic signihcant flow 
rate reduction is observed, partly explaining the “increased” or “effective” viscosity reported 
in the literature. Other effects, like anisotropy, will also play a role in the change in effective 
transport properties. For polar molecular systems like water /c ~ 3 — 4 nm, and the coupling 
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must be considered on these length scales. For the non-polar fluids studied here Ic is below 
1 nm and the coupling effect is very small in most situations. 

(ii) In general any fluid response can be described phenomenologically through a trans¬ 
port kernel that incorporates the spatial and temporal correlation effects. A method for 
calculating the shear viscosity kernel was presented. This showed that for non-highly vis¬ 
cous fluids the Newtonian limit is reached on length scales of a few nanometers. This is 
in agreement with the Bocquet-Charlaix criterion if the decay time for the stress autocor¬ 
relation function is applied as the relaxation parameter. Importantly, non-local effects are 
not present in simple flows where the strain-rate is linear with respect to spatial coordinate, 
which is the case for Couette and Poiseuille flows. For non-linear flows the non-local response 
significantly affects the fluid stress for strain-rate variations on the atomic length scale. 

(hi) For highly confined fluids molecular alignment phenomena and molecular deformation 
can occur along with the fluid layering, see Fig. 8. Simple classical continuum theory does 
not include or account for such complex fluid structure. It would be interesting to investigate 
this is more detail, for example, using the theory for liquid crystals [70, 81]. 

In conclusion, continuum theory is applicable even on the nanoscale if the relevant phys¬ 
ical processes are modelled appropriately. 
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SUPPORTING INFORMATION 


BALANCE EQUATIONS AND THE EXTENDED NAVIER-STOKES EQUA¬ 
TION 

A macroscopic field variable A(r, t) is be given by the corresponding microscopic quantity 
ttiit) associated with molecule i through [1] 

i 

where (5(r) is the Dirac delta function and rj(t) is the molecular center-of-mass position. In 
our treatment A(r, t) can be a scalar quantity {e.g. the mass density) or a vector quantity 
{e.g. the momentum density). In Fourier space we have for A{r,t) 

p P i*00 _ 

J'[A{r,t)] = A(k,t) = JJJ A(r, dr = ^ . (60) 

From now on we do not write the time and position dependencies of the free variables (right 
hand side of the equations) explicitly unless it provides important information. The rate of 
change in Fourier space is given by 

^I(k, = Y1 + ai{-ik ■ Vi) \ (61a) 

= (^1 - ■ Ti - ^(k ■ Tif .. .^ , (61b) 

where Vj is the center-of-mass velocity. The following identity has been applied, where a is 
a scalar or a vector 

a(b ■ c) = (b ■ c)a = b ■ (ca) . ( 62 ) 

In the situation where a is a vector the product ca is the dyadic between two vectors, c and 
a, and the resultant is the second rank tensor with components {ca)ij = Cittj, i,j = x,y,z. 
The dyadic is not commutative, i.e., ca ^ ac, unless ca is symmetric. The dyadic is also 
called the outer product and sometimes written as c 0 a. If c and a are parallel the dyadic 
is symmetric. To hrst order in wave vector Eq. (61b) becomes 

^A(k,t) = ^(1 - fk ■ r*)^ - fk ■ ^ ViOi (small k). (63) 
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The particle velocity can be decomposed into a thermal (or ’’peculiar”) term c* and advective 
term u(rj,t) by 

Vi(t) = Cj(t) + u(ri,t) . (64) 

u(rj,t) is the mass averaged fluid velocity in a region around r^. This region is equivalent 
to the fluid particle volume encountered in the traditional continuum description [2] as was 
briefly discussed in the introduction. We assume that the thermal and advective velocities 
are uncorrelated and the thermal motion is, by deflnition, conserved with = 0 , 

where rrii is the mass. From this Eq. (63) becomes 

d ~ 

—A{k,t)='H[ai] (small k), (65) 

where "H is a linear operator given by 

'^N = ■ r*)^ - ik ■ - ik ■ ^u{ri,t)ai . (66) 

i i i 

H is henceforth refered to as the “microscopic hydrodynamic operator” (MH-operator) since 
it describes the microscopic interpretation of the fluid’s dynamics in the hydrodynamic 
regime of small wave vectors. The first term describes the rate of change of the microscopic 
quantity in question, the second and third terms are then the thermal and advective contri¬ 
butions to the dynamics, respectively. We proceed and apply Eqs. (65) and (66) to express 
the balance equations of mass, linear momentum and angular momentum. 


A. Mass balance 


The mass density held p(r, t) is defined directly from Eq. (59) as [3] 

p{r,t) = - r^) . 

i 

The MH-operator acting on rrij gives 

'H[mi]ikA) = ^(1 - ik ■ r*)"^ “ ■ '^rriiCi - ik ■ ^miu(ri,t) 

i i i 

= -ik ■ 


(67) 


(68a) 

(68b) 


since the mass of each particle is constant and rrijCj = 0 . The dynamics in Fourier space 
is thus simply 


A 


■p(k,t) 


—ik ■ 


rriiUi 


i 


(small k). 


(69) 
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Relating to the general balance equation, Eq.(l) in the manuscript, in the case of zero 
production term, one has 0 = 1, and no surface forces 

=-V-pu . (70) 

In Fourier space this yields 

^p(k,t) =-ik-^ (71) 

using that J^[—V ■ (pu)](k, t) = —ik ■ pu, which follows from partial integration of Eq.(60). 
In the limit of small wave vector we thus identify Eq. (69) as the microscopic interpretation 
of the continuous mass balance equation. 

B. Linear momentum 

In general, there is no microscopic dehnition of the streaming velocity u in the form 
of Eq. (59) [3]. Rather, the linear fluid motion is given through the momentum density 
J(r, t) = p(r, t)u(r, t). Note, Hansen and McDonald [1] dehne a current from the microscopic 
velocities and Eq. (59), which is then the correct mass weighted averaged streaming velocity 
for single species fluids. If is the linear momentum of particle i we have 

J(r, t) = p(r, t)u(r, t) = ^ rniVi6{r - r^) . (72) 

i 

Applying the MH-operator we have 

'H[miVi]i^^,t) = ^(1 - ik ■ Vi)mi^ - ik ■ - ik ■ ^ miu(ri, t)vj (73a) 

i i i 

= ^(1 - ik ■ ri)Ff* - ik ■ ^ ruiCiCi - ik ■ ^ miu(ri, t)u(rj, t) (73b) 

i i i 

as it can be shown [4] that the cross terms miu(rj, t)cj = miCju(ri, t) = 0 since the 
thermal and advective velocities are uncorrelated. 

F*"* is the total force acting on molecule i. This force is decomposed into the force on 
i due to interactions with all other molecules denoted henceforth by Fj and external forces 
F|^*, i.e., F(°* = Ff^* + Fj. The first term on the right hand side of Eq. (73b) is 

5^(1 - ik . r.)F“ = 5^(1 - ik. r.)Fr' - *k • r.F. (74) 
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using Newton’s third law X^jFj = 0 and the identity Eq. (62). Substituting this into Eq. 
(73b) and rearranging one arrives at 

(small k), (75) 

i i / 

where 

5„ = 5^(l-ik.r.)Fr‘. (76) 

i 

We recognize Eq. (75) as the balance equation for linear momentum in Fourier space. Im¬ 
portantly, the last term on the right hand side gives one possible microscopic interpretation 
of the pressure tensor P in the zero wave vector limit [4]. Let V be the system volume. Then 
according to Eq. (1) in the manuscript and (75) we have for the zero wave vector pressure 
tensor 

VP(t) = P(k = 0, t) = ^ rUiCiCi + ^ rjFj. 

i i 

This is the Irving-Kirkwood interpretation [5] and will be discussed in the following. 

The conhgurational part of the pressure tensor is the term X^jTjFj and can be written 
in a different form assuming pairwise additive interactions only. Using Newton’s third law 
once again, Fij = —Fji, where Fj^ is the force on i due to interactions with j, we get 

TjFj = Fij = J'pFjj , (77) 

i i j i j>i 

where = r* — r^. The Irving-Kirkwood pressure tensor then reads 

rp(i) = E TTliCiCi -|- EE 

i i j>i 

At this point we emphasize the difference between the atomic and molecular formalisms. 
If particles i and j represent point-masses (including atoms in the molecules), rather than 
structured molecules, the vectors Vij and Fij are parallel, the dyadic rj^F^ is symmetric. 
This implies that the pressure tensor is symmetric. In the case of structured molecules, 
however, the force on molecule i due to j is given by 

F*j(^) = Fiaj/3 , (79) 

a&i fiej 

where indices ia and j/9 represent atom a in molecule i and atom /? in molecule j, see Fig. 
11. The force F^ needs not be parallel to the vector connecting the molecules’ center-of- 
masses and the conhgurational part of the pressure tensor is not generally symmetric. There 


d 

—pu(k,t)-Fik-^miu(ri,t)u(ri,t) = a-^-ik- 
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FIG. 11. Schematic illustration of the interactions between two molecules i and j. The atoms are 
represented by open circles and chemical bonds between atoms by straight lines. Filled circles are 
the center-of-masses. The force acting on i due to j, Fjj, needs not to be parallel with the vector 
Tij connecting the two center-of-masses. 

is no ambiguity between the two formalisms. They are trivially equivalent for point-mass 
particles, and for structured molecules one can be expressed in terms of the other and the 
mass dispersion tensor [4, 6]. The two formalisms provide different information and we use 
the molecular one here. We use the term ’atom’ loosely: it represents a point-mass spherical 
particle, and can be a single atom or a united group of atoms like the methyl group. 

In the molecular formalism the pressure tensor decomposes into a sum of the symmetric 

s a s a 

P and anti-symmetric P parts such that P =P -|- P with 

P(i) = i(P + Pq and P (i) = i(P - PP . (80) 

Following Evans and Morriss [3] we here use the stack notation to indicate the symmetric 
and anti-symmetric parts. The symmetric part can be further decomposed into a sum of 

OS 

the equilibrium pressure peq, the viscous pressure If, and the trace-less symmetric tensor P 

[3] 

S OS 

P (t) = {Peq + n)I+ p , (81) 

where I is the second rank identity tensor and peq + If = Tr(P)/3. The anti-symmetric part 
of the pressure tensor has zero diagonal components and three independent off-diagonal 

ad a a a 

components. It is therefore often written as a (pseudo) vector dual P= {Pyz, Pzx, Pxy) of 

a 

the tensor. Here the subscripts indicate the tensor components. From the dehnition of P, 
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Eq. (80), one has for P 

Ojd » ^ 

i 

= ^ Fi X Fi = ^ ^ r^ X Fij . (82) 

i i j>i 

Here i, j and k are the unit vectors along the x-, y- and x-axes, respectively. From Eq. (82) 
it is clear that the anti-symmetric pressure is associated with the torque on i about j and 


gives rise to a change of orbital angular momentum, he., that of the center-of-mass motion 
of the molecules. Importantly, when the pressure tensor has a non-zero anti-symmetric 
part, the orbital angular momentum is a non-conserved quantity [7]; we address this in the 
following. Conservation of orbital angular momentum is sometimes used as an argument 
for a symmetric pressure, see for example Aris [8], but note that later Aris also allows for 
the possibility of asymmetry. Finally, notice that the kinetic part of the pressure tensor 
is a symmetric dyadic and does not enter the anti-symmetric part. 


C. Orbital and spin angular momenta 

Consider now the case when no external forces are present. Following the approach for 
the linear momentum we write the orbital angular momentum density p(r, t)L(r, t) as 

p(r, t)L(r, = ^ ^ Pi)<^(r “ r*) , (83) 

i i 

where L(r, t) is the orbital angular momentum per unit mass, and L* = r* x p* is the orbital 
angular momentum of molecule i. Note, the angular momentum is dehned with respect to 
some specihc choice of coordinate system. The MH-operator acting on Lj gives 

H[Li](k,t) = ^(1 -ik ■ r,)-^ - fk ■ ^ CiLi - ik ■ ^ u(rj, f)Li (84a) 

i i i 

= N - ik ■ ^ r^Ni - ik ■ ^ c^Li - ik ■ ^ u(ri, t)Li, (84b) 

ill 

where N* = r* x Fj is the torque on molecule i and N is the sum of torques 

a ^ ad 

N = ^riXFi = 2FP . (85) 

i 

The torque N does not sum to zero as the orbital angular momentum is a non-conserved 
quantity [7]. Rather it is the total angular momentum, i.e., the sum of orbital and spin (or 
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intrinsic) angular momenta, which is conserved. For point-masses spin angular momentum 
is not present (or meaningful) as there is no moment of inertia. In this case the orbital 

a 

angular momentum is conserved and the torque N must be zero, that is, P= 0. The rate of 
change of orbital angular momentum is then given by 

Q ___ 

—pL(k,f)-h fk ■ ^u(ri,t)Li = N - fk ■ ^ (cjLj-h r^Ni) (small k). (86) 

i i 

This is the balance equation for the orbital angular momentum in Fourier space. 

The spin angular momentum density [9] p(r,t)S(r, t) is given by 

p(r,t)S(r,t) = ^ Si(5(r - r*) = ^ | x j 5(r - r*). (87) 

i i \ ia / 

S(r,t) is the spin angular momentum per unit mass, Si(t) = ^ P*a is the spin 

angular momentum of molecule i, and Rja is the vector from the center-of-mass to atom a, 
see Fig. 11. Applying the MH-operator and following the procedure in Eq. (84) above, the 
rate of change is to lowest order in wave vector 

^pS(k, t)-h fk ■ ^ UjSj = M - ik ■ ^ (ciSi + TiMj) (small k), (88) 

i i 

where Mj = X^j^Ria x Fjo is the sum of torques on i about the center-of-mass and M = 
Mj. Equation (88) is the balance equation for the spin angular momentum in Fourier 
space. 

The last term in Eq. (88) dehnes the zero wave vector Irving-Kirkwood couple tensor 
[10] VQ(t) = Q(k = 0,t). The flux of spin angular momentum associated with the couple 
tensor can be important for flows in extreme conhnements. After volume averaging Q can 
be written as 

VQ{t) = 

i i j>i 

where 

^ Rio X ^ 

FjajP (90) 

id jl3 

is the torque, specihcally the couple force, on molecule i due to interaction with all atoms 
P in molecule j, see Fig. 11. The couple tensor is not symmetric and can be decomposed 
into symmetric and anti-symmetric parts just as above for the pressure tensor, that is. 


05 a 

Q = Qi-\- Q -l- Q 


(91) 
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ad a a a 

Also, one has a vector dual Q= {Qyz)Qzx^Qxy) for the anti-symmetric part of the couple 
tensor. 

The total angular momentum is a conserved quantity. This means that in the limit of 
zero wave vector 

+ /^) = 0, z.e., M + N = 0 (92) 

at 

and so we have from Eqs. (82) and (85) 

E ad 

r, X F, = -2F P . (93) 

i 

Above we discussed the mass and momentum balance equations in the framework of 
the microscopic picture, which is easiest done in Fourier space. We will now return to the 
corresponding real space formulations, but write the balance equations in a slightly different 
form compared to Eq. (1) in the main manuscript. In real space Eq. (70) is the mass 
balance equation. Using the product rule on the right hand side and re-arranging 

^ = -p(V . u) , (94) 


where, in general, the material derivative D/Dt is dehned by 

_ = _ + (V^), 


(95) 


Likewise, from Eqs. (75) and (88) we get the relevant momentum and spin balance equations 


in real space after application of Eq. (94) 

Du _ 05 ad 

p— = cTu - V ■ {{peg + n)I-h P) -h Vx P (96a) 

DS ad os ad 

P-— = o-s - 2 p -V ■ (gi+ Q) + Vx Q , (96b) 


where erg represents the external production term for the spin angular momentum. In Eqs. 

a ad 

(96) we have applied the identity V- A= — Vx A (A = P or Q). Interestingly, the term 

ad 

—2 P in Eq. (96b) can be regarded as a production term even in the absence of an external 
torque. Now, Eqs. (94), (96a) and (96b) give the relevant balance equations in the limit of 
small wave vectors or, equivalently, large length scales. 

Also, from Eqs. (96) one immediately sees that by ignoring the anti-symmetric part 
of the pressure tensor we obtain the classical momentum balance equation for the linear 
momentum. As discussed above this applies to point-mass structure-less fluids. Also, one 
observes that the dynamics of the spin angular momentum is then determined by the couple 
forces alone and the spin is a conserved quantity. 
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D. The extended Navier-Stokes equations 


The treatment above is general. We shall now focus on systems where the molecules can 
be well approximated as uni-axial rigid body molecules. Uni-axial molecules are defined 
as molecules where the principal moment of inertia tensor has two nonzero components, 
denoted Ip, the third one being zero as the rotation around the main molecular axis is not 
associated with any inertia. This includes di-atomic and linear molecules such as carbon 
dioxide, but other molecules are also well approximated as uni-axial and rigid, for example, 
butane as we will see below. The treatment in this section is not as detailed as above, and 
we refer the reader to Refs. [7, 11-13]. 

It is convenient to describe the spin dynamics in the principal frame of reference where 
the moment of inertia tensor is constant. From the Euler equation of motion for rigid bodies 
[14] one can show that in this frame and for uniaxial molecules the left hand side of Eq. 
(96b) becomes 


DS Dfl 
P-RI = P^- 


Dt 


Dt 


(97) 


where I = 2Ip/3 [11, 12] and f2 is the angular velocity. Note, this does not hold in general 
where nonlinear coupling between the angular velocity components are present. The general 
situation, however, can be included into the theory. 

For isotropic systems the following linear local constitutive relations are applied, see Refs. 
7 and 10, 


n = -? 7 ^(V ■ u) 


OS 

p 


ad 

p 


os / 2 

-2?7o(Vu) = -? 7 o f (Vu -h uV) - -Tr(Vu)I 
—riri^ XU — 217) 


Q = -C.(V-17) 

os os f 2 

Q = -2Co(V17) = -Co i (V12 + 17V) - -Tr(V17)I 

ad 

Q = -Cr(V X 17) . 


(98a) 

(98b) 

(98c) 

(98d) 

(98e) 

(98f) 


Here rjo, and are the shear, bulk and rotational viscosities, respectively, and Co; Cv and Cr 
are the corresponding spin viscosities. Usually, a thermodynamic force is defined through the 
gradients of the held variables, however, in Eq. (98c) we recognize the thermodynamic force 
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V X u — 217 as twice the difference between the classical prediction of the angular velocity, 
= iV X u, and the actual one. This force is denoted the sprain rate [10]. Substituting 
Eqs. (98) into Eqs. (96) and using Eq. (97) we arrive at the extended Navier-Stokes (ENS) 
equations [10] 

= CTJ - Vpeq + {r]v + po/3 - ?7r)V(V ■ u) + (?7o + r]r)V^u + 2T]rV X Q (99a) 

D17 n 

= cTg + 2?7r(V Xu — 217) + ((C^ + Co/3 — Cr)V(V ■ 17) + (Co + Cr)V 17 (99b) 

5 

using the divergence rule V ■ (Vu) = + |V(V ■ u). 

We now present values for the transport coefficients entering Eqs. (98) for the molecular 
fluids studied here. Methane is modelled as a point-mass molecule where the coupling is 
irrelevant, he., we let = 0. Butane and water are not uni-axial. For the butane model 
used the principal moment of inertia components are Ix = 2.2 x 10“^° m^, ly = 1.8 x 10“^'^ 
m^ and Iz = 0.26 x 10“^° m^, that is, the molecule has one major axis and two almost 
identical minor axes and we can expect the theory to hold reasonably well. In Ref. 15 the 
water molecule was considered as a dipolar rotator with an effective moment of inertia of 
I = 8.4 X 10“^^ m^. This approach will be adopted here. In Table I the relevant state points 
and transport coefficients are listed. In the treatment below the dynamics is decomposed 
into transverse shear and longitudinal bulk modes. For linear momentum the focus is on the 
shear mode, and the bulk viscosity is therefore not listed in the table. The coefficients are 
evaluated from independent equilibrium MD simulations as prescribed in Refs. 10, 15-17. 
The values for the dumbbell model are given in reduced MD units of length scale a, energy 
scale e, and mass m, see Ref. 16. 


Molecule 

p [kg m ^] 

T[K] 

po [mPa-s] 

rjr [mPa-s] 

Co + Cr [kg m s"^] 

I [m2] 

Methane 

460 

164.4 K 

0.27* 

- 

- 

- 

Dumbbell** 

0.4477 

4.0 

0.60 

0.083 

0.22 

1/6 

Butane 

582.3 

288 K 

0.14 

0.013 

4.0 X 10-24 

1.3 X10-20 

Water 

996.3 

298.7 

0.7 

0.17 

2.1 X 10-21 

8.4 X 10-22 


TABLE I. * From Ref. 18. ** Quantities given in reduce MD units. 

State points, transport coefficients and moments of inertia for the systems studied. 
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